Question: Does a continuous wetland to upland gradient exist? Or, can a catena-type model be derived using wetlands and uplands in a data-driven approach?

hydropedology hillslope model
hydropedology hillslope model

Using Hubbard Brook as an example



hbdem <- rast("UplandWetlandGradient/data/dem1m.tif")
plot(hbdem)

hb <- vect("UplandWetlandGradient/data/hbef_wsheds/hbef_wsheds.shp")

lines(hb)

nwi <- vect("UplandWetlandGradient/data/HU8_01070001_Watershed/HU8_01070001_Wetlands.shp") |>
    project("EPSG:26919")
hbnwi <- crop(nwi, ext(hbdem))
names(hbnwi)
## [1] "ATTRIBUTE"  "WETLAND_TY" "ACRES"      "Shape_Leng" "Shape_Area"
plot(hbdem)
plot(hbnwi, "WETLAND_TY", type = "classes", add = T)

hbgeo <- vect("UplandWetlandGradient/data/hbef_bedrock/hbef_bedrock.shp")

plot(hbgeo, type = "classes")


# 1 meter LiDAR-derived, Filtered Hydro-enforced Digital
# Elevation Model (DEM)
hb1mlpns <- rast("UplandWetlandGradient/data/hydem1mlpns.tif")


plot(hb1mlpns)

Get training data - wetland and upland points

randupl <- spatSample(hbdem, 2000, as.points = TRUE, xy = TRUE)

hbnwi_buff <- buffer(hbnwi, 50)

hbupl_pts <- terra::mask(randupl, hbnwi_buff, inverse = T) |>
    mutate(class = "UPL") |>
    select(-dem1m)


plot(hbupl_pts, type = "points")
lines(hbnwi_buff)

hbwet_pts <- spatSample(hbnwi, size = 500) |>
    mutate(class = "WET") |>
    select(class)

plot(hbwet_pts, type = "classes")


hbpts_all <- rbind(hbupl_pts, hbwet_pts)
hbpts_ext <- terra::extract(hbdem, hbpts_all, bind = T)
hbpts <- hbpts_ext |>
    dplyr::filter(!is.na(dem1m))
writeVector(hbpts, "UplandWetlandGradient/data/derived_data/hbpts_raw.gpkg",
    overwrite = TRUE)

plot(hbdem)
plot(hbpts, "class", type = "classes", add = T)

Make terrain metrics

hbslp_3 <- SlpAsp(hbdem, w = c(3, 3), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_3.tif",
    overwrite = T)

hbslp_27 <- SlpAsp(hbdem, w = c(27, 27), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_27.tif",
    overwrite = T)

hbslp_81 <- SlpAsp(hbdem, w = c(81, 81), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_81.tif",
    overwrite = T)

hbtpi_3 <- TPI(hbdem, w = c(3, 3), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_3.tif",
    overwrite = T)

hbtpi_27 <- TPI(hbdem, w = c(27, 27), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_27.tif",
    overwrite = T)

hbtpi_81 <- TPI(hbdem, w = c(81, 81), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_81.tif",
    overwrite = T)

hbcurv_3 <- Qfit(hbdem, w = c(3, 3), metrics = c("meanc", "profc",
    "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_3.tif",
    overwrite = T)

hbcurv_27 <- Qfit(hbdem, w = c(27, 27), metrics = c("meanc",
    "profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_27.tif",
    overwrite = T)

hbcurv_81 <- Qfit(hbdem, w = c(81, 81), metrics = c("meanc",
    "profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_81.tif",
    overwrite = T)

Slope at different scales

plot(hbslp_3)
plot(hbslp_27)
plot(hbslp_81)

TPI at different scales

plot(hbtpi_3)
plot(hbtpi_27)
plot(hbtpi_81)

Curvature at different scales

plot(hbcurv_3)
plot(hbcurv_27)
plot(hbcurv_81)

super_stack <- c(hbslp_3, hbslp_27, hbslp_81, hbcurv_3, hbcurv_27,
    hbcurv_81, hbtpi_3, hbtpi_27, hbtpi_81)

hbpts_extract <- terra::extract(super_stack, hbpts, bind = TRUE,
    xy = TRUE)
names(hbpts_extract) <- (c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81", "x", "y"))
hbpts_extract$class <- as.factor(hbpts_extract$class)

writeVector(hbpts_extract, "UplandWetlandGradient/data/derived_data/hbpts_extract.gpkg",
    overwrite = TRUE)
hist(hbpts_extract$slp_3)
hist(hbpts_extract$slp_27)
hist(hbpts_extract$slp_81)

Slope at 3m vs. 81m

plot(hbpts_extract$slp_3, hbpts_extract$slp_81)
abline(0, 1, col = "red")

TPI at 3m vs. 81m

plot(hbpts_extract$tpi_3, hbpts_extract$tpi_81)
abline(0, 1, col = "red")

Mean Curvature at 3m vs. 81m

plot(hbpts_extract$meancurv_3, hbpts_extract$meancurv_81)
abline(0, 1, col = "red")

Set up random forest model

library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:spatialEco':
## 
##     combine
set.seed(11)

hbpts_exdf <- as.data.frame(hbpts_extract)

# Validation Set
train.index <- as.vector(sample(c(1:nrow(hbpts_exdf)), 0.7 *
    nrow(hbpts_exdf), replace = F))

train <- hbpts_exdf[train.index, c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81")]

test <- hbpts_exdf[-train.index, c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81")]
library(caret)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
## Loading required package: lattice
rf_model <- randomForest((class) ~ ., ntree = 500, na.action = na.omit,
    importance = TRUE, data = train)
plot(rf_model)

rf.test <- predict(rf_model, newdata = test)

caret::confusionMatrix(rf.test, test$class)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction UPL WET
##        UPL 497  16
##        WET   9 102
##                                           
##                Accuracy : 0.9599          
##                  95% CI : (0.9414, 0.9739)
##     No Information Rate : 0.8109          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8663          
##                                           
##  Mcnemar's Test P-Value : 0.2301          
##                                           
##             Sensitivity : 0.9822          
##             Specificity : 0.8644          
##          Pos Pred Value : 0.9688          
##          Neg Pred Value : 0.9189          
##              Prevalence : 0.8109          
##          Detection Rate : 0.7965          
##    Detection Prevalence : 0.8221          
##       Balanced Accuracy : 0.9233          
##                                           
##        'Positive' Class : UPL             
## 
vip::vip(rf_model)

AUC/ROC

Predict maps

pred_stack <- c(hbdem, super_stack)


names(pred_stack) <- (c("dem", "slp_3", "slp_27", "slp_81", "meancurv_3",
    "prof_curv_3", "plan_curv_3", "meancurv_27", "prof_curv_27",
    "plan_curv_27", "meancurv_81", "prof_curv_81", "plan_curv_81",
    "tpi_3", "tpi_27", "tpi_81"))
hbwip <- predict(pred_stack, rf_model, type = "prob", filename = "UplandWetlandGradient/data/derived_data/hbwip.tif",
    overwrite = TRUE)

plot(hbwip$WET)
lines(hb)

PDF/kernel density plot

hbwip_trainpred <- predict(rf_model, newdata = train, type = "prob")[,
    2]
hist(hbwip_trainpred, 50)

plot(density(hbwip_trainpred, na.rm = T, from = 0, to = 1))

hb_pedons <- vect(read.csv("UplandWetlandGradient/data/HBEF_pedon_locations.csv"),
    geom = c("easting", "northing"), crs = "EPSG:26919")

hb_pedons_ext <- terra::extract(hbwip, hb_pedons, bind = T, method = "bilinear")
hb_pedons_ext$hpu <- factor(hb_pedons_ext$hpu, ordered = T)

plot(hbwip$WET)
plot(hb_pedons, "hpu", type = "classes", add = T)

ggplot(hb_pedons_ext, aes(x = reorder(hpu, desc(WET)), y = WET)) +
    geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), adjust = 1.5) +
    geom_point(aes(color = hpu), alpha = 0.5, position = position_jitter(seed = 11))
## Warning: Removed 4 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).

ggplot(hb_pedons_ext, aes(x = reorder(pm, desc(WET)), y = WET)) +
    geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), adjust = 1.5) +
    geom_point(aes(color = pm), alpha = 0.5, position = position_jitter(seed = 11))
## Warning: Removed 4 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Removed 4 rows containing missing values (`geom_point()`).

plot(hbwip$WET)
plot(hb_pedons_ext[hb_pedons_ext$WET >= 0.5, ], "hpu", type = "classes",
    add = T)

w3_pedons <- vect(read.csv("UplandWetlandGradient/data/hbef_w3_lat_pedons/W3-lateral-weathering-pedons.csv"),
    geom = c("easting", "northing"), crs = "EPSG:26919")

w3_pedons_ext <- terra::extract(hbwip, w3_pedons, bind = T, method = "bilinear")

w3 <- hb[hb$WS == "WS3"]
w3_wip <- crop(hbwip, w3, mask = T)
plot(w3_wip$WET, legend = F)
plot(w3_pedons, "hpu", type = "classes", add = T, legend = T)
lines(hbgeo)
lines(hbnwi)

ggplot(w3_pedons_ext, aes(x = hpu, y = WET)) + geom_violin(draw_quantiles = c(0.25,
    0.5, 0.75), adjust = 1.5) + geom_point(aes(color = hpu),
    alpha = 0.5, position = position_jitter(seed = 11))
## Warning: Removed 1 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

ggplot(w3_pedons_ext, aes(x = water_table, y = WET)) + geom_point(aes(color = hpu),
    alpha = 0.5, position = position_jitter(seed = 11))
## Warning: Removed 30 rows containing missing values (`geom_point()`).